# Cargar tmap
library(tmap)
# Colección de paquetes de Tidyverse
library(tidyverse)
# Estilos para ggplot2
library(ggthemes)
# Paletas de colores de RColorBrewer
library(RColorBrewer)
# Paletas de colores de viridis
library(viridisLite)
# Gráficos interactivos
library(plotly)
# Manejo de datos vectoriales
library(sf)
# Manejo de datos raster
library(terra)
# Manejo de datos raster
library(raster)
# Mapas interactivos
library(leaflet)
# Acceso a datos en GBIF
library(rgbif)
# Datos geoespaciales
library(geodata)
# Modelado de distribución de especies
library(dismo)
# Manejo datos Maxent
library(rJava)
# Cargar la librería de widgets
library(htmlwidgets)Aves Guanacaste
Cargar librerías
Contexto
Los datos de este estudio se obtuvieron de tres fuentes distintas:
eBird: que es una plataforma de ciencia ciudadana que recopila, almacena y analiza datos de aves de personas observadoras en todo el mundo.
WorldClim: proporciona datos climáticos a nivel mundial, estos incluyen variable climáticas como tempertura, precipitación, radiación solar, humedad relativa y otros factores importantes para la modelación climatológica y ecológica.
HydroSHED: que tiene datos acerca del flujo del agua, acumulación de flujo, cuencas hidrológicas, río y lagos y permiten el análisis detallado sobre los sistemas hidrográficos a escala global.
Definimos el área de Guanacaste como área de estudio debido a que es un territorio que cuenta con 2 sitios Ramsar: el Embalse Arenal (8,317 ha) y el humedal del PN Palo Verde los cuales son sitios de reproducción y alimentación para una gran cantidad de especies de aves acuáticas, migratorias y residentes. Estos humedales funcionan como refugio de especies en peligro de extinción y constituye una de las zonas de anidamiento más grandes del país (SINAC,2020 y 2024)
Nuestro objeto de estudio son las aves: Dendrocygna autumnalis, Spatula discors, Jabiru mycteria, Mycteria americana, Eudocimus albus. Estas especies tienen una vida vinculada a los humedales, pues sus poblaciones y ciclos de vida dependen de estos ecosistemas (alimentación, reproducción, anidación, etc). Todas son residentes a excepción Spatula discors que tiene un periodo de migración al norte, entre septiembre y abril, para reproducirse (Barchiesi et al., 2021).
Cargar datos aves y área de estudio
# Lectura de un archivo TXT con registros de presencia de aves en Costa Rica
puntos <- read.delim(
"C:/Users/laura/OneDrive - Universidad de Costa Rica/Monitoreo_ambiental/Aves_Palo_Verde/Ebird/Ebird data Guanacaste/ebd_CR-G_200701_202412_unv_smp_relSep-2024/ebd_CR-G_200701_202412_unv_smp_relSep-2024.txt"
)
# Cargar shape Guanacaste
guanacaste <- st_read("C:/Users/laura/OneDrive - Universidad de Costa Rica/Monitoreo_ambiental/Monitoreo_Palo_Verde/Guanacaste/delimitshapes/guanacaste.shp",
quiet = TRUE)Filtrar especies seleccionadas
# Filtrar las especies de aves de interés
especies_seleccionadas <- puntos |>
filter(SCIENTIFIC.NAME == "Dendrocygna autumnalis" |
SCIENTIFIC.NAME == "Spatula discors" |
SCIENTIFIC.NAME == "Jabiru mycteria" |
SCIENTIFIC.NAME == "Mycteria americana" |
SCIENTIFIC.NAME == "Eudocimus albus")Definición espacial: sistema referenciado de coordenadas.
# Convertir a objeto espacial usando columnas de coordenadas
especies_seleccionadas <- st_as_sf(
especies_seleccionadas,
coords = c("LONGITUDE", "LATITUDE"), # Nombres exactos de las columnas
remove = FALSE, # conservar las columnas de las coordenadas
crs = 4326 # Sistema de referencia WGS84 (EPSG:4326)
)
# Asignación de un CRS al objeto guanacaste
guanacaste <-
guanacaste |>
st_transform(4326)#Dar formato de fecha a la columna correspondiente
especies_seleccionadas$OBSERVATION.DATE <- as.Date(especies_seleccionadas$OBSERVATION.DATE, format = "%Y-%m-%d")Graficar registros de presencia por año por especie
Podemos observar la cantidad de registros de las especies dependen del interés de las y los observadores, por lo que especies interesantes o difíciles de conseguir presentan muchos registros, como el Mycteria americana, a diferencia de aquellas que son más especies más comunes y por tanto más fáciles de encontrar como Dendrocygna autumnalis
# Gráfico ggplot2
grafico_ggplot2 <-
especies_seleccionadas |>
st_drop_geometry() |>
group_by(year = year(OBSERVATION.DATE), SCIENTIFIC.NAME) |>
summarize(n = n()) |>
ggplot(aes(x = year, y = n, color=SCIENTIFIC.NAME)) +
geom_line() +
geom_point(
aes(
text = paste0(
"Año: ", year, "\n",
"Cantidad de registros: ", n
)
)
) +
ggtitle("Cantidad de registros de presencia por año") +
xlab("Año") +
ylab("Cantidad de registros de presencia") +
labs(caption = "Fuente: eBird") +
theme_minimal()
# Gráfico plotly
ggplotly(grafico_ggplot2, tooltip = "text") |>
config(locale = 'es')Gráfico de barras anual
Se muestra la cantidad de individuos por especie registrados cada año. Se puede observar que los patos presentan mayor número de individuos (Dendrocygna autumnalis y Spatula discors).
# Agrupo y cuento los individuos por año
df <- especies_seleccionadas |>
group_by(SCIENTIFIC.NAME, year = year(OBSERVATION.DATE)) |>
summarise(
individuos = sum(as.numeric(OBSERVATION.COUNT), na.rm = TRUE)
) # Elimino nans
# Gráfico de barras
graf_barras <- df |>
ggplot(aes(x = as.factor(year), y = individuos, fill = as.factor(SCIENTIFIC.NAME))) +
geom_bar(
stat = "identity",
aes(
text = paste0(
" Individuos por año: ", round(individuos, 2)
)
)
) +
ggtitle("Número de Individuos por Año") +
xlab("Año") +
ylab("Cantidad de Individuos") +
labs(caption = "Fuente: Datos de eBird", fill = "Especies") +
theme_minimal()
# scale_fill_manual(values = colores, name = "Año") # Pongo los colores definidos y nombre a la leyenda
# Gráfico interactivo con Plotly
ggplotly(graf_barras, tooltip = "text") |>
config(locale = 'es')Cargar los archivos de Clima
# Definir el directorio de los datos climáticos
ruta_clima <- "C:/Users/laura/OneDrive - Universidad de Costa Rica/Monitoreo_ambiental/Aves_Palo_Verde/Modelos/climate/bios/"
# Listar los archivos .tif en el directorio
archivos_clima <- list.files(ruta_clima, pattern = "\\.tif$", full.names = TRUE)
# Cargar todos los archivos raster
bios <- rast(archivos_clima)
names(bios) [1] "wc2.1_30s_bio_1" "wc2.1_30s_bio_10" "wc2.1_30s_bio_11" "wc2.1_30s_bio_12"
[5] "wc2.1_30s_bio_13" "wc2.1_30s_bio_14" "wc2.1_30s_bio_15" "wc2.1_30s_bio_16"
[9] "wc2.1_30s_bio_17" "wc2.1_30s_bio_18" "wc2.1_30s_bio_19" "wc2.1_30s_bio_2"
[13] "wc2.1_30s_bio_3" "wc2.1_30s_bio_4" "wc2.1_30s_bio_5" "wc2.1_30s_bio_6"
[17] "wc2.1_30s_bio_7" "wc2.1_30s_bio_8" "wc2.1_30s_bio_9" Recortar área de estudio
# Definir la extensión del área de estudio
area_estudio <- ext(
min(especies_seleccionadas$LONGITUDE) - 5,
max(especies_seleccionadas$LONGITUDE) + 5,
min(especies_seleccionadas$LATITUDE) - 5,
max(especies_seleccionadas$LATITUDE) + 5
)
# Recortar las variables bioclimáticas al área de estudio
clima <- crop(bios, area_estudio)Cargar datos de humedales
Nosotras utilizamos la capa de Global Lakes and Wetlands Database: GLWD_v2_delta_area_ha_x10, que indica el porcentaje de cobertura de humedales, que para éste análisis es adecuado puesto que las aves y sus hábitos de vida están relacionados a ellos.
ruta_humedal <- "C:/Users/laura/OneDrive - Universidad de Costa Rica/Monitoreo_ambiental/Aves_Palo_Verde/GLWD_v2_delta_combined_classes_tif/GLWD_v2_delta_combined_classes/GLWD_v2_delta_area_ha_x10.tif"
archivo_humedales <- rast(ruta_humedal)
# Recortar y resamplear 'humedales' para que tenga la misma extensión y resolución que 'clima'
humedales <- archivo_humedales |>
crop(ext(clima)) |>
resample(clima)Unir datos de clima y humedales en un mismo ráster
# Unir las capas de clima y de humedales
datos <- c(clima, humedales)Entrenamiento del modelo MaxEnt
Se crea un DataFrame con las coordenadas de Longitud y Latitud
# Crear dataframe con columnas de longitud y latitud
coordenadas_presencia <- data.frame(
decimalLongitude = especies_seleccionadas$LONGITUDE,
decimalLatitude = especies_seleccionadas$LATITUDE
)
# Eliminar coordenadas duplicadas
coordenadas_presencia <- unique(coordenadas_presencia)Crear semilla para garantizar aleatoriedad reproducible
# Establecer una "semilla" para garantizar que la selección aleatoria sea reproducible
set.seed(123)
# Cantidad de registros de presencia
n_presencia <- nrow(coordenadas_presencia)
# Con sample(), se selecciona aleatoriamente una proporción (ej. 0.7)
# de los índices de los datos de presencia para el conjunto de entrenamiento
indices_entrenamiento <- sample(
1:n_presencia,
size = round(0.7 * n_presencia)
)
# Crear el subconjunto de entrenamiento utilizando los índices seleccionados
entrenamiento <- coordenadas_presencia[indices_entrenamiento, ]
# Crear el subconjunto de evaluación con los datos restantes
evaluacion <- coordenadas_presencia[-indices_entrenamiento, ]Cambiar formato a ráster y ejecutar el modelo
# Los datos de clima deben convertirse al formato que usa el paquete raster
# debido a que es este el que acepta el paquete dismo
datos <- raster::stack(datos)
# Ejecutar el modelo
modelo_maxent <- maxent(x = datos, p = entrenamiento)
# Aplicar el modelo entrenado a las variables climáticas
# para generar un mapa de idoneidad del hábitat
prediccion <- predict(modelo_maxent, datos)Evaluacion del modelo MaxEnt
# terra::extract() extrae los valores del raster de predicción
# en las coordenadas de evaluación
# eval_pres almacena los valores de idoneidad predichos
# en los puntos de evaluación de presencia
eval_pres <- terra::extract(
prediccion,
evaluacion[, c('decimalLongitude', 'decimalLatitude')]
)
# Generar puntos aleatorios dentro del área de estudio definida.
# Estos puntos se asumen como ausencias de la especie.
ausencias <- randomPoints(mask = datos, n = 1000)
# eval_aus almacena los valores de idoneidad predichos
# en los puntos de ausencia
eval_aus <- terra::extract(
prediccion,
ausencias
)
# Generar estadísticas de evaluación del modelo
resultado_evaluacion <- evaluate(p = eval_pres, a = eval_aus)Curva ROC
La curva ROC y los valores del área bajo la curva (AUC) dan un buen ajuste, puesto que el modelo está prediciendo de buena manera la distribución de las aves que ingresamos. El valor de AUC corresponde a 0.99 y en esta ocasión y en un intento anterior 0.993, lo cual es muy cercamp a 1 por lo que el modelo encaja de buena manera con la distribución real.
# Datos para graficar la curva ROC
datos_roc <- data.frame(
FPR = resultado_evaluacion@FPR,
TPR = resultado_evaluacion@TPR,
Umbral = resultado_evaluacion@t
)
# Valor AUC
auc <- resultado_evaluacion@auc
# Gráfico ggplot2
grafico_ggplot2 <-
ggplot(
datos_roc,
aes(
x = FPR,
y = TPR,
u = Umbral
)
) +
geom_line(
color = "blue",
size = 1
) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "red") +
labs(title = paste("Curva ROC (AUC =", round(auc, 3), ")"),
x = "Tasa de falsos positivos (FPR)",
y = "Tasa de verdaderos positivos (TPR)") +
theme_minimal()
# Gráfico plotly
ggplotly(grafico_ggplot2) |>
config(locale = 'es')Graficar mapa de distribución continua
En este mapa se muestra en un rango de 0 a 1, donde cero es ausencia y presencia, la posiblidad de que las aves estén presentes. Esto da un mapa con degradado donde encontramos áreas con alta probabilidad de presencia cercanas a los humedales y zonas altas y baja probabilidad en zonas bajas y con menos precipitación.
# Paleta de colores del modelo
colores_modelo <- colorNumeric(
palette = c("white", "black"),
values(prediccion),
na.color = "transparent"
)
# Paleta de colores de temperatura
colores_temperatura <- colorNumeric(
# palette = "inferno",
# palette = "magma",
palette = rev(brewer.pal(11, "RdYlBu")),
values(datos$wc2.1_30s_bio_1),
na.color = "transparent"
)
# Paleta de colores de precipitación
colores_precipitacion <- colorNumeric(
# palette = "viridis",
# palette = "YlGnBu",
palette = "Blues",
values(datos$wc2.1_30s_bio_12),
na.color = "transparent"
)
# Paleta de colores de precipitación
colores_humedales <- colorNumeric(
# palette = "viridis",
# palette = "YlGnBu",
palette = "YlGnBu",
values(datos$Band_1),
na.color = "transparent"
)
# Crear una paleta de colores categórica para las especies
especies_unicas <- unique(especies_seleccionadas$SCIENTIFIC.NAME)
paleta_especies <- colorFactor(
palette = brewer.pal(length(especies_unicas), "Set1"),
domain = especies_unicas
)
# Mapa
leaflet() |>
addTiles(group = "Mapa general") |>
addProviderTiles(
providers$Esri.WorldImagery,
group = "Imágenes satelitales"
) |>
addProviderTiles(
providers$CartoDB.Positron,
group = "Mapa blanco"
) |>
addRasterImage( # capa raster de temperatura
datos$wc2.1_30s_bio_1,
colors = colores_temperatura, # paleta de colores
opacity = 0.6,
group = "Temperatura",
) |>
addRasterImage( # capa raster de precipitación
datos$wc2.1_30s_bio_12,
colors = colores_precipitacion, # paleta de colores
opacity = 0.6,
group = "Precipitación",
) |>
addRasterImage( # capa raster de humedales
datos$Band_1,
colors = colores_humedales, # paleta de colores
opacity = 0.6,
group = "Humedales",
) |>
addRasterImage( # capa raster del modelo de distribución
prediccion,
colors = colores_modelo,
opacity = 0.6,
group = "Modelo de distribución",
) |>
addCircleMarkers(
# capa de registros de presencia (puntos)
data = especies_seleccionadas,
stroke = F,
radius = 3,
fillColor = ~paleta_especies(SCIENTIFIC.NAME),
fillOpacity = 1,
popup = paste(
paste0("<strong>Nombre Científico: </strong>", especies_seleccionadas$SCIENTIFIC.NAME),
paste0("<strong>Cantidad de individuos: </strong>", especies_seleccionadas$OBSERVATION.COUNT),
paste0("<strong>Fecha de avistamiento: </strong>", especies_seleccionadas$OBSERVATION.DATE),
paste0("<strong>Fecha de avistamiento: </strong>", especies_seleccionadas$COMMON.NAME),
paste0("<a href='", especies_seleccionadas$GLOBAL.UNIQUE.IDENTIFIER, "'>Más información</a>"),
sep = '<br/>'
),
group = "Registros de aves"
) |>
addLegend(
title = "Temperatura",
values = values(datos$wc2.1_30s_bio_1),
pal = colores_temperatura,
position = "bottomleft",
group = "Temperatura"
) |>
addLegend(
title = "Precipitación",
values = values(datos$wc2.1_30s_bio_12),
pal = colores_precipitacion,
position = "bottomleft",
group = "Precipitación"
) |>
addLegend(
title = "Humedales",
values = values(datos$Band_1),
pal = colores_humedales,
position = "bottomleft",
group = "Precipitación"
) |>
addLegend(
title = "Modelo de distribución",
values = values(prediccion),
pal = colores_modelo,
position = "bottomright",
group = "Modelo de distribución"
) |>
addLayersControl(
# control de capas
baseGroups = c("Mapa general", "Imágenes satelitales", "Mapa blanco"),
overlayGroups = c(
"Temperatura",
"Precipitación",
"Humedales",
"Modelo de distribución",
"Registros de aves"
)
) |>
hideGroup("Temperatura") |>
hideGroup("Precipitación") |>
hideGroup("Humedales")